package src.montecarlo;

import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.IOException;
import java.util.Random;

import src.percolation.Ej1;



public class Ising2D {

	//parameters of the system
	private double magneticField;
	private double couplingFactor;
	private double kT;

	//variables of the system
	private boolean[][] grid;
	private double[] probabilityTransitions;
	private double[] energyTransitions;
	private int large;
	private double energy;
	private double magnetization;
	private int spinsDifferential;
	private int spinPairs;
	private Random generator;
	
	public Ising2D(int large, double magneticField, double couplingEnergy,double temperature) {
		
		if(large<0 || magneticField<0 || couplingEnergy<0 || temperature<0){
			System.out.println("System parameters must be positive.");
			return;
		}
		
		this.large = large;
		this.magneticField = magneticField;
		this.couplingFactor = couplingEnergy;
		this.kT = temperature;
		
		this.grid = new boolean[large][large];

		generator = new Random();
		magnetization = populateGrid();
		energy = getSystemEnergy();
		loadTransitionProbabilitiesAndEnergies();
	}
	/**
	 * Useful: this is needed only once, as there is a very limited set of possible on state flips
	 * @param couplingEnergy
	 * @param temperature
	 * @return
	 */
	private void loadTransitionProbabilitiesAndEnergies() {

		//duplicated because of magnetic degeneration 
		probabilityTransitions = new double[10];
		energyTransitions = new double[10];
		
		
		
		double[] coupledEnergy = new double[5];
		
		coupledEnergy[0] = 8 * couplingFactor;
		coupledEnergy[1] = 4 * couplingFactor;
		coupledEnergy[2] = 0 * couplingFactor;
		coupledEnergy[3] =-4 * couplingFactor;
		coupledEnergy[4] =-8 * couplingFactor;
		
		double[] magneticEnergy = new double[2]; 
		
		magneticEnergy[0] = 2* magneticField;
		magneticEnergy[1] =-2* magneticField;
		
		for(int i=0;i<2;i++){
			for(int j=0;j<5;j++){
				energyTransitions[i*5+j] = (coupledEnergy[j]+magneticEnergy[i])/(large*large);
				probabilityTransitions[i*5+j] = Math.exp(-(energyTransitions[i*5+j]/kT));
			}
		}
	}

	public Ising2D loadConfiguration(File config){
		
		int large = 0;
		double magneticField = 0;
		double couplingEnergy = 0;
		double temperature = 0;
		//also populate grid accordingly
		
		Ising2D ising = new Ising2D(large,magneticField,couplingEnergy,temperature);
		//aca hay que popular la red como esta en la configuracion
		
		
		return ising;
	}
	
	
	public void saveConfiguration(String config) throws IOException{
		
		FileOutputStream fos = new FileOutputStream(new File(config));
		//header and variables
		fos.write(("/large\tmagnetic field\tcoupling energy\ttemperature\tenergy\tmagnetization\n").getBytes());
		fos.write((large+"\t"+magneticField+"\t"+couplingFactor+"\t"+kT+"\t"+energy+"\t"+magnetization+"\n").getBytes());
		
		//complete state
		for(int i=0;i<large;i++){
			for(int j=0;j<large;j++){
				if(grid[i][j])	fos.write("1".getBytes());
				else			fos.write("0".getBytes());
			}
			fos.write("\n".getBytes());
		}
	}
	
	/**
	 * Only useful the first time.
	 * Later on, we should take care whether magnetization is updated.
	 */
	private double getSystemEnergy() {
		
		spinPairs = 0;
		
		for(int i=0;i<large;i++){
			for(int j=0;j<large;j++){
				
				if(j<large-1){
					if(grid[i][j]==grid[i][j+1]) spinPairs++;
					else spinPairs--;					
				}else if (j==large-1){//periodic boundary conditions
					if(grid[i][j]=grid[i][0]) spinPairs++;
					else spinPairs--;
				}
				
				if(i<large-1){
					if(grid[i][j]==grid[i+1][j]) spinPairs++;
					else spinPairs--;
				}else if(i==large-1){//periodic boundary conditions
					if(grid[i][j]==grid[0][j]) spinPairs++;
					else spinPairs--;
				}
			}
		}
		return (-couplingFactor * spinPairs - magnetization*spinsDifferential)/(large*large);
	}

	/*
	 *Calculates the magnetization of the system while populating the grid 
	 */
	private double populateGrid() {
		
		spinsDifferential=0;
		
		for(int i=0;i<large;i++){
			for(int j=0;j<large;j++){
				boolean state = generator.nextBoolean();
				grid[i][j] = state;
				if(state)
					spinsDifferential++;
				else
					spinsDifferential--;
			}
		}
		return spinsDifferential*magneticField/(large*large);	
	}

	/*
	 * Update system energy and magnetization
	 */
	private boolean flipOneSpin(int x, int y){
		
		boolean state = grid[x][y];
		//energy update:
		int nettedNeighboursSpin = 0;
		//left neighbor
		if(x>0){
			if(grid[x-1][y])		nettedNeighboursSpin++;
			else					nettedNeighboursSpin--;
		}else if(x==0){
			if	(grid[large-1][y])  nettedNeighboursSpin++;
			else 				    nettedNeighboursSpin--;	
		}
		//right neighbor	
		if(x<large-1){
			if(grid[x+1][y]) 		nettedNeighboursSpin++;
			else		   		    nettedNeighboursSpin--;
		}else if(x==large-1){
			if(grid[0][y])	 		nettedNeighboursSpin++;
			else			 		nettedNeighboursSpin--;
		}
		//down neighbor
		if(y>0){
			if(grid[x][y-1])		nettedNeighboursSpin++;
			else					nettedNeighboursSpin--;
		}else if(y==0){
			if(grid[x][large-1]) 	nettedNeighboursSpin++;
			else			 	 	nettedNeighboursSpin--;
		}
		//up neighbor
		if(y<large-1){
			if(grid[x][y+1]) 		nettedNeighboursSpin++;
			else			 		nettedNeighboursSpin--;
		}else if(y==large-1){
			if(grid[x][0])	 		nettedNeighboursSpin++;
			else			 		nettedNeighboursSpin--;
		}

		double probability = getMockProbability(state, nettedNeighboursSpin);
				
		//Metropolis-MonteCarlo acceptance procedure
		double threshold = generator.nextDouble();

		if(probability>threshold){//accept the new state
			acceptNewState(x,y,state,nettedNeighboursSpin);//E & M updated inside
			return true;
		}
		return false;
	}
	
	public void shakeTheSystemOnce(){
		
		int rejectedFlips = 0;
		
		//we will flip L^2 spins randomly, and update energy and magnetization accordingly on each step.
		for(int i=0;i<large*large;i++){
			int x = Math.abs(generator.nextInt()%large);
			int y = Math.abs(generator.nextInt()%large);
			if(!flipOneSpin(x, y))//update E and M inside
				rejectedFlips++;
		}	
	}
	public void shakeTheSystemManyTimes(int n){
			for(int k=0;k<n;k++) shakeTheSystemOnce();
	}

	public double getEnergy() {
		return energy;
	}

	public double getMagnetization() {
		return magnetization;
	}
	
	private void acceptNewState(int x, int y, boolean state,int deltaSpinPairs){
		
		if(grid[x][y]){//it was up, now down
			spinsDifferential--;
			magnetization -= magneticField*2/(large*large);
		}
		else{
			spinsDifferential++;
			magnetization += magneticField*2/(large*large);
		}
		
		energy += getDeltaEnergy(state,deltaSpinPairs);
		spinPairs += deltaSpinPairs;
		
		//finally, invert it
		grid[x][y] = !grid[x][y];
	}
	
	private double getDeltaEnergy(boolean state, int nettedNeighbors) {
		
		int i;
		int j=0;
		
		//this part of the energy is owed to magnetic field
		if(state) i=1;
		else i=0;
		
		//this part of the energy because of coupling spins
		switch(nettedNeighbors){
		case 4:
			if(state)	j=0;
			else;		j=4;
			break;
		case 2:
			if(state)	j=1;
			else		j=3;		
			break;
		case 0:
			j=2;
			break;
		case -2:
			if(state)	j=3;
			else		j=1;
			break;
		case -4:
			if(state)	j=4;
			else		j=0;
			break;
		}
		return energyTransitions[i*5+j];
	}
	/**
	 * This is a mock probability as it may return a value >1.
	 * It's supposed to be used against a double random generator [-1,1]
	 * and accept a new state in the case of having a bigger probability
	 * proba>1 equals deltaEnergy<0
	 * 
	 * @param state
	 * @param nettedNeighbors
	 * @return
	 */
	private double getMockProbability(boolean state,int nettedNeighbors){
		
		int i;
		int j=0;
		
		if(state) i=1;
		else i=0;
		
		switch(nettedNeighbors){
		case 4:
			if(state)	j=0;
			else;		j=4;
			break;
		case 2:
			if(state)	j=1;
			else		j=3;		
			break;
		case 0:
			j=2;
			break;
		case -2:
			if(state)	j=3;
			else		j=1;
			break;
		case -4:
			if(state)	j=4;
			else		j=0;
			break;
		}
		return probabilityTransitions[i*5+j];
	}
}
